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Abstract 

The generalization of the simplest equation method to look for exact 
solutions of systems of nonlinear differential equations is presented. The 
exact solutions of NDE systems describing the evolution of two interacting 
populations in two cases (both populations have the low critical density 
or low critical density is typical for only one of populations) are obtained. 

Introduction. 

Systems of nonlinear differential equations describe a lot of processes in 
different fields of science. Since XIX century it became clear that nonlinear 
phenomena are no less important than the linear ones. Therefore, the theory of 
nonlinear differential equations began its development. 

Usually we use the approximate methods to solve the problems with the 
systems of nonlinear differential equations, but exact solutions are very useful 
too. Exact analytical solutions allow to determine the features of solutions 
behavior under some initial and boundary conditions and verify the numerical 
calculations. 

In [1] the simplest equation method to look for exact solutions of nonlinear 
differential equations was presented (see also [2-5]). It summarizes a number of 
methods that were developed earlier, including the tanh-function method [6-9], 
the Jacobi elliptic function method [10,11] and the Weierstrass function method 
[12,13]. 

Here we generalize the simplest equation method for systems of ordinary 
differential equations. This method applies two main concepts: 1) the idea 
of the simplest nonlinear differential equation that has lesser order than the 
equations of the system studied and 2) accounting of possible singularities of 
the solutions of system studied. 

The outline of this paper is as follows. The generalization of the simplest 
equation method to look for exact solutions of systems of nonlinear differential 
equations is given in section 1. Sections 2 and 3 are devoted to the exact 
solutions of NDE systems describing the evolution of two interacting populations 
in two cases 1) if both populations have the low critical density (section 2) and 
2) if low critical density is typical for only one of populations (section 3) [14]. 



1 The simplest equation method to look for ex- 
act solutions of systems of nonlinear differen- 
tial equations. 

Consider the system of m nonlinear differential equations in the polynomial 
form 

Mi[yi,y2, ■ ■ -ym] ^ 0, i=l,2,...n (1) 
and the simplest equation 

E[Y] = (2) 

with order lesser than the ones of equations Mi. The solution of simplest equa- 
tion is assumed to be known. 

For example, we can use as the simplest equation the Riccati equation, the 
equation for the Jacobi elliptic function, the equation for the Weierstrass elliptic 
function, and so on. 

The unknown functions t/j, j — 1, 2, ... m are supposed to be expressed via 
the solution of the simplest equation Y. 

Let us describe the simplest equation method to look for exact solutions of 
systems of nonlinear differential equations. It contains three basic steps. 

At the first step we determine the singularity orders of unknown functions 
?/j , j = 1, 2, . . . m by the analysis of the dominant terms of system ^ (this step 
coincides with the first step of the Painleve test). 

At the second step the solutions of ODE system ([T]) are represented as se- 
lected simplest equation solution expansions, usually in the polynomial form 

yj^F,{Y)^Y.^MY) (3) 

k 

The terms used in expansions are selected in consideration of singularity orders 
of simplest equation ^ and singularity orders of the studied system solutions, 
that were determined at the first step. 

At the third step we substitute expansion ^ into the system studied. After 
consideration of used simplest equation properties, coefhcients of various Y- 
function powers are set equal to zero. So we obtain the algebraic system, which 
solution determine coefficients Ajk and, if it is necessary, the restrictions to the 
parameters of initial system. 

Further we consider the applications of the simplest equation method to 
the systems of nonlinear differential equations that are used in ecology at the 
description of interacting species. 
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2 Exact solutions of the nonlinear differential 
equations system describing the evolution of 
two interacting species with low critical den- 
sities. 

Evolution of two interaction populations with low critical densities can be 
described by the reaction-diffusion system in the form 

du ^ (Pu 

— — Ui-r^ + aiu[Ki — u)(u — Li) + El — eiuv 

at dx^ ^^-j 

dv d^v 

— = -^2^^ + a2V[K2 ~ V)[V - L2) + E2- S2UV 

Here u = u{x, t) and v = v{x, t) are the densities of populations, t means time, 
X is spatial variable. 

Type of interactions between species is determined by signs of parameters Si 
(competition at ei, £2 > 0, symbiosis at ei, £2 < 0, "predator-prey" or "parasite- 
host" at £i£2 < 0). 

Parameters Di, a;, Ki, Li and Ei are determined for each population sep- 
arately and describe the evolution of population without species interactions 
{i = 1, 2). Coefficient Di characterizes the velocity of population diffusion by 
influence of random motion [15], define the growth rate, Ki is the environ- 
ment capacity, Li is the low critical density (if population is under minimum 
viable level it become extinct), Ei describe the external influence on ecosys- 
tem [14,16,17]. 

System of equations ^ can be reduced to the Burgers-Huxley equation 
[18-20] by linear substitution if and only if the diffusion velocities of both pop- 
ulations coincide, i.e. Di — D2, and two additional conditions on parameters of 
system (|4]) hold. However, it is important to have exact solutions for different 
diffusion coefhcients. 

Substituting the variables 



2Co^ , / 2Co^D2 Co 



= \ F^"' z=—{x-CQt) 



and denoting 



Pi = 



haiDi aiDiKiLi , aiDi^ 

-^{Li+Ki), p2 = — —2 , El ^^ly^^^ 



2a2Di ,r ^r.^ a2DiK2L2 ^, ^ a2Di' 



, Di , _ 2D2 , _ I 2Di- 

a -7^, £i-£iW — —J, £2-^21 



n ' 1 '-\ 2 ' '-2 ^ 2 n2 

L>2 y ^260 V Q^i^o D2 

omitting strokes we obtain the dimensionless form of system ^ in traveling 
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wave variables 





du 


dz^ ' 


dz 


d^v 


^dv 


dz"^ ' 


dz 



(5) 



2v + qiv — q2V + E2 — e2uv = 

Use the simplest equation method. 

The dominant terms of system ([S]) are the following 



(6) 



This truncated system has solution u{z) = si/z, v{z) = S2/Z, where si and S2 
can take on values ±1 independently. Therefore, solutions of system ([5]) have 
first order singularities. 

As the simplest equation we use the Riccati equation in the form 

^ = -y2 + p^Y + Po (7) 
dz 

where Y = Y{z) and Po, Pi are constants. 

The Riccati equation belongs to the class of exactly solvable ones, its solu- 
tions have first order singularities. Nonhomogeneous solutions of equation ^ 
can be written as 

• At P2 + 4Po > 



Y{z) = \{Pi + \/P? + 4Po tanh(^p2 ^ 4p^(^ _ J (g) 



Yiz) = ^ ( ^1 + V^i + "^^0 coth(^p2 + 4Po(z - zo)) ) (9) 



At P2 + 4Po = 



Yiz) = ^ + (10) 
2 z — zq 



• At + 4Po < 

y{^) = ^ (^1 - \/-(^' + 4^") tan(y^-(p2 + 4Po)(z - zq))) (H) 

Y{z) ^l[Pi + ]/-iP? + ^Po) C0t(^-(P2 + 4Po)(z - zo))^ (12) 

Here Zq is the constant of integration. Note, that function ^ is the only 
one bounded on complete number scale, however other solutions can be used on 
bounded intervals. 

As solutions of system ^ arc to have first order singularities, we use the 
substitution 

_ ^ ^ Y 

u{z) ^ AqY + Ai + A2y 

v{z) = BoY + Bi + B2Y 
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Taking into the account the properties of function Y ( i.e. the Riccati equation 
([7])), the equivalent form of this substitution is 

u{z) = AoY + Ai + ^ 

BP ^13) 

v{z) = BoY + B, + 

Nonhomogeneity condition for sohition (|13p is 

Al + AlP^ ^ 0, El + BlP^ ^ (14) 

Statement. 

We can assume Aq ^ without loss of generality. 

Proof. Suppose that Aq — Q. If A2P0 = then the nonhomogeneity condi- 
tion mi) fails, therefore Pq 7^ 0, ^2 7^ 0. Denote Y Pa/Y. 

% = {P^m. = -y^Y, = -y^i-Y' + PiY + Po) = 

= -^-PiY+Po = -Y^ + PiY + Po (15) 

so function Y satisfy the Riccati equation in the form ([7]) with Pi = —Pi. 
Therefore solution at = can be presented as 

u{z) = AqY + Ai 

v{z) = BoY + Bi + ^ 
Y 

where Ao ~ A2, Bo = B2 and B2 — Bo- □ 
Substituting (|13p in system ([5]), taking into acctount 



Y, = + p^Y + Po 



- _ sp^y^ + (p2 „ 2Po)Y + PoPi 



(16) 



and setting coefficients of various y-function powers equal to zero, we get the 
system of fourteen algebraic equations. Consider its solutions. 

Equations generated by the largest and the least powers of Y are 

Ao{Al-l)^0, A2{Al-l) = 0, Po(Bo-l)-0, 52(^2 - 1) = 

so we can sort possible solutions by values of parameters Ao, A2, Bo and B2. 

• Case A2= B2= 0, A^ = B^ = 1. 

Here substitution p3|) is invariant under linear transformations, so we 
can consider Pi = without loss of generality. Therefore we define the 
parameters of solution 

Aopi - BoEi - 1 
Al — — 



Bi 



&Ao 
BoQi - Aoe2 



GBo 

Po = -l/12si{-e2+Pi)BoAo + l/12ei (d- l)Po 
+ 1/I2pi2 - l/12eigi - 1/12 - l/2p2 



(17) 
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and coefficients of initial system, that allow this solution 

92 - 1/6 qi^ - 1/6 pi^ + l/Qeiqi +P2 + 1/6 - l/6£2Pi - 1/6 - 
- 1/6 £2 {d - 1) Ao - 1/6 £i {d - 1) Bo + (1/6 £iPi - 1/6 qi£2) B^Aq 
El = Ai (2 Ai^ - piAi + Bi£i + Pa) - ^o^o 
E2 = Bi (2 Bi^ - qiBi + e2Ai + 52) - dBoPo 



• Case ^2 = Bo = 0, = S| = 1. 

Here system ([5|) has solution in the form (|13p with 



at condition that 



qi = B2il + d - piAo) 
q2 = (l + d + 3p2-pi (d+l)Ao)/3 
^1 = Po^o (£1^2 -Pi -I) 
E2 = P0B2 (£2^0 + Pi-d) 



(18) 



Ai = Bi=Q 
Pi = [piAo - l)/3 (19) 
Po = {piAa-2+pi^ - %P2)/IS 



(20) 



• Case A2 = 0, Al = Bl = Bl = 1. 
Here 

Ai = 
Pi = (piAo - 1 - £lBo)/3 
Pi = (gi+B2(3Pi-d))/6 
Po = (A+Pi'-P2-eiPi)/2 

£2 = Ao(d(PoB2 - 1) - 3Pi(l + B0B2)) 

£;i = AoPo(£iP2 - Bi - 1) 
£;2 = Po ((-2 qi + £2^oBo + 12 Pi) P0P2 - (Po + B2) - Pi (Po - P2)) - 

+ 2Pi3 + 92Pi -91P1' 
92 = Pi^ - 6P1 + 2Pigi - dPi - 2Po(l + 3P0P2) 
qi = P2(rf- 3Pi) + 12PidAoPo/£2 

(21) 

• Case P2 = 0, = = P§ = 1. 
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Figure 1 



Here 

Bi = 
Pi - (giBo-d-e2^o)/3 
Ai-(pi+^2(3Pi-l))/6 
Po = {dPi+Pi^ -q2~e2Ai)/2 
ei = So(AoA2 - 1 - 3Pi(l + A0A2)) 
El = Po {{-2pi + e2AoBo + 12 Ai) A0A2 - {Aq + A2) - Pi {A^ - A2)) + 

+ 2Ai^ +P2Ai-piAi^ 
E2 = BoPo{e2A2 -Pi-d) 

P2 = P2 _ g^2 ^ 2AiPi - Pi - 2Po(l + 3Ao^2) 

pi = ^2(1 - 3Pi) + UPidAoBo/ei 

(22) 

Other solutions of algebraic system either don't satisfy the conditions of 
positiveness of pi, qi and d or reduce to solutions sited above. 

If we consider solutions bounded on complete number scale only, the so- 
lutions obtained describe two types of populations behavior: 1) simultaneous 
monotone change of both populations densities (fig. [T]) and 2) change of one of 
populations state by act of the interactions with the solitary wave of the other 
population (fig. [2]) 

Note, that behavior of the second type is generated by solutions in case 3 
at B2 = —-Bo and in case 4 at A2 = — Ao. Thus solution given on fig. [5] is 
obtained at pi — 4.8, p2 = 2.8, d — 0.8, Si — 3.1, qi = 7.7, q2 = 2.45, £2 = —1.6, 
El = 0.234, £^2 = and has the form 
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Figure 2 



23 V V 
u{z) = ^ + tanh ( -^(z - ^o) 



V6T 


tanh 1 


V6T 


20 




20 


%/6T 




%/6T 


tanh 


20 




20 



, , 23 V61 ^ / V61, .\ 117 
z;(z) = — — — tanh \z — zq 



20 20 \ 20 ^ 5(23 + V6Ttanh(# 



where zq is the arbitrary constant (at fig. [2] it assumed to be zero). It corre- 
sponds to the situation when the predators rise results in the prey extinction, 
but predators can subsist in the bounded domain only. 



Exact solutions of the nonlinear differential 
equations system describing the evolution of 
two interacting species, one of which is char- 
acterized by low critical density. 

Consider the system of nonlinear differential equations in the form 



du a'u . , , , 

— ~ Ui-r^ + 2aiu(Ai — u)[u — Li) + hi — eiuv 
at ax'' 



(23) 
av a V 

— = '^'^ a2v[K2 - V) + E2- e2uv 

In contrast to system (01), it describes the situation, when only one of popula- 
tions is characterized by low critical density, and the other one in the absence 
of interacting species and external influence complies with the Fcrhulst-Pearl 
equation [14,16,17,21]. The meanings of parameters of system coincide with 
the ones declared in section 2. 
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Substituting the variables 



u = 2{Ki + Li)u, V =2 — —[Ki+Lij v, 



a2Di ' ' ' V 2ai {Ki + Li 

and denoting 



2(Xi + Li)2' 4ai(Li + is:i)3' Diaa' 

^1 ^ vT^ ;~Tj ^2 — £2- 



V2aTi>^(A'l aiD2{Ki+Li) 

omitting strokes we obtain the dimensionless form of system (|23p in travehng 
waves variables 



cPu ^ du ^ ^ ^ ^ 

-r^ +Ci- + u- pu + El - siuv = 

az'' dz 

—i; + dCi + qv + E2- e2uv = 

az'' dz 

The dominant terms of system are the following 

d^u „ o 

— -77 — 2u — eiuv = 



(24) 



(25) 



The solution of this truncated system is u{z) ~ si-^/l — 3ei/z, w(z) = 6/2^, 
where si can take on value ±1. 

So the solution of system (f24|) we are looking in the form 



u{z) = AoY + Ai + 



v{z) = SqJ-^ + BiY + B2 



^ 2 (26) 



Y r2 

where y = Y{z) is the solution of the Riccati equation ([7]) (here we take into 
account that function u has first order singularities and function v has the second 
order singularities). 

The nonhomogeneity condition for solution in the form ()26p is the following 

+ A2^Po^ ^ 0, Bo^ + Bi^Po^ ^ (27) 

Substituting transformation ((26|) in system (|24|) . taking into account the 
Riccati equation ([7]) and setting coefficients of various y-function powers equal 
to zero, we obtain sixteen algebraic equations. Equations generated by the 
largest and the least powers of Y are 

Ao(2-2V-£iSo) = 
^2^0^(2 - 2A2^ - £1^4) = 
So(Bo-6) = 
BiPo\Bi - 6) = 
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We can sort three cases by values of parameters Aq, A2, Bq and B4, which 
specify aU possible solutions of system in the form ((25)) under natural 
conditions of parameters of initial system. Consider them sequentially. 

• Case Aq^ = 1 - 3ei, Bq = 6, = S4 0. 

As a result we have that B3 — Q. At A2 ^ B4 — transformation 
is invariant under linear transformations, so we assume Pi = without 
loss of generality. Then we can find sequentially determine the solution 
parameters 



Ai 



Pi = -3(£2Ao + 2Cid)/5 
(1 - 3£i) (5 + 3£i£2) - AqCi (5 -6eid) 
30(1 - 2£i) 



£2^ - 3 £i£2^ + 125 £2^1 - 125 g + 17 Ci^d^ 
200(12C7id + £2A„) + (2^) 

, „ 8Ci2d2 + 27£i£22-9£22 



200(12Cid + £2Ao) 
P2 = (6 g - Pi^ - CirfPi - 48 Po - 6 £2^1 - £2AoPi)/12 

and coefficients of the initial system, at which this solution exists 
p = 2 Ai - 6 Ai^ _ 2 Po - £iP2 



^0 
+ 2 

P2 = £2AiP2 - 12 Po^ - CidBiPo + B2^ - qB; 



Pi = pAx - Ai^ - CiAoPq + eiAiB2 + 2 Ai^ (^^^ 



Case = 1 - 3£i, Al = 1, Pq = 6, P4 = 0. 

Here the solution parameters and coefhcients of system (12^ are to satisfy 
the relations 

Pi = Al (2 Ai^ - Al +p) - Po (Pi + Ci) Ao - 

^ Po (-eiPi + Ci - Pi) A2 + 2 Po (6 yli - 1) A2AQ 
E2 = e2^2PoSi - 12 Po^ - PiPiPo - CidPiPo 
p = -6A0A2P0-2P0-C1P1 -6Ai2 + Pi2 + 2Ai 
g = -4 Pi^ + 1/2 Pi Pi + £2^1 + 1/6 Pi^ + 1/6 CidPi + 
+ 1/6 £2AoPi - 2 CidPi + 8 Po 
Pi (-6 CidPi - 18 Pi^ + 3 PiPi + Pi^ + CidPi + £2AoPi) 



Pi 



36(Pi + 6 Pi - £2^2 + 2 Cid) 
Pi = -3/5£2Ao -6/5Cid-6Pi 
Al = 1/6+ (1/2 Pi - l/6Ci)A2 
P2 = P3 = 

Ci (6 £id - 5) Ao - 5 Ci (-1 + 2 £1) A2 - £1 (9£i£2 + 5 - 3 £2) 



15(l-2£i)(Ao + A2) 

(30) 
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Figure 3 



and condition 

(CiPi-6Ai2-p+2Ai+Pi2-2Po) Ao- 

-6 (l-3ei)A2Po-ei(^iSi+6A2Po) =0 (31) 

(which can be reduced to the cubic equation in C\) must hold. 

• Case Aq^ = 1 - 3ei, A2 = -Aq, Bo = B4 = 6. 
Solution with parameters 

2d -5 



Ai 



6(d-5)' 



R _ 9^V2 , 5d(2d-5) 
' 200 ' ^ 16(d-5)2' 



o _ 2d 6d 



Pi = 
is allowed at 



3(d-5) 



^° = 400^^ ~ 



19d(2d- 5) 
288(d-5)2 



d{4d - 25) 



400 



5(4d-F3)(2d-5) 
96(d-5)2 ' 

3d(2d- 5) 



3^1^ 2 

20 ' ^ 8(d-5)2 



£1 



Eo = 



_5_ 
6d 
351d^ 



£2 



20d 



3(d-5)' 



^1=0, 



40000 



4 21d3(2d-5)^2 9d2(2d-5)2 
Ci — 



64(d - 5)2 



256(d - 5)4 



(32) 



(33) 



Solutions, given in this section, can demonstrate more complicated behavior 
in comparison with ones obtained in section 2. Some examples are presented 
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Figure 4 



at figures [21 HI and [5l Thus at fig. 2] density of one of populations under the 
influence of the other population diffusion falls before the equilibrium state 
mounts. 

Conclusion. 

In this paper the generalization of simplest equation method to look for exact 
solutions of nonlinear differential equations is described. 

Using the simplest equation method, we obtain the exact solutions of systems 
describing the evolution of two interacting species in two cases 1) when both 
populations are characterized by low critical densities and 2) when low critical 
density is important for one of populations only and the other one follows the 
Ferhulst-Pearl law. 

All the solutions obtained describe the spatial transitions between the steady 
states of the system studied. 
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